home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QFLTI.C < prev    next >
C/C++ Source or Header  |  1996-04-02  |  36KB  |  2,224 lines
  1. /* ------------------------------------------------------------- *
  2.  *                            qflti.c  *
  3.  *            QFLOAT                     *
  4.  *                                 *
  5.  *      EXTENDED PRECISION FLOATING POINT ROUTINES         *
  6.  *                                 *
  7.  * Prototype        Purpose                  *
  8.  * ---------        -------                  *
  9.  * asctoq(string, q)    ascii string to q type             *
  10.  * dtoq(&d, q)        DEC double precision to q type         *
  11.  * etoq(&d, q)        IEEE double precision to q type      *
  12.  * e113toq(&d, q)    128-bit long double to q type         *
  13.  * e64toq( &e, q )    80-bit IEEE long double to q type     *
  14.  * ltoq(&l, q)        long integer to q type             *
  15.  * qabs(q)        absolute value                 *
  16.  * qadd(a, b, c)    c = b + a                 *
  17.  * qclear(q)        q = 0                     *
  18.  * qcmp(a, b)        compare a to b                 *
  19.  * qdiv(a, b, c)    c = b / a                 *
  20.  * qifrac(x, &l, frac)    x to integer part l and q type fraction  *
  21.  * qldexp(x, n)     multiply x by 2^n             *
  22.  * qinfin(x)        set x to infinity, leave its sign alone  *
  23.  * qmov(a, b)        b = a                     *
  24.  * qmul(a, b, c)    c = b * a                 *
  25.  * qmuli(a, b, c)    c = b * a, a has 16 significant bits     *
  26.  * qisneg(q)        returns sign of q             *
  27.  * qneg(q)        q = -q                     *
  28.  * qnrmlz(q)        adjust exponent and mantissa         *
  29.  * qsub(a, b, c)    c = b - a                 *
  30.  * qtoasc(a, s, n)    q to ASCII string, n digits after decimal*
  31.  * qtod(q, &d)        convert q type to DEC double precision     *
  32.  * qtoe(q, &d)        convert q type to IEEE double precision  *
  33.  * qtoe113(q, &d)    convert q type to 128-bit long double     *
  34.  * qtoe64( q, &e )    q type to 80-bit IEEE long double     *
  35.  *                                 *
  36.  *                                 *
  37.  *      EXTENDED PRECISION MATHEMATICAL FUNCTIONS         *
  38.  *                                 *
  39.  *                                 *
  40.  * qexp( x, y )     y = exp( x )                 *
  41.  * qfloor( x, y )    y = largest integer not greater than x     *
  42.  * qlog( x, y )     y = log( x ) [natural logarithm]     *
  43.  * qpow( x, y, z )    z = x^y [x raised to the y power]     *
  44.  * qrand( q )        q = pseudorandom number in [0,1)     *
  45.  * qremain( a, b, c )    c = remainder after dividing b by a.     *
  46.  * qround( x, y )    y = nearest integer to x         *
  47.  * qsqrt( x, y )    y = sqrt( x )                 *
  48.  * qsrand( u )        initialize seed of qrand w/ unsigned int *
  49.  * qtanh( x, y )    y = tanh( x ) [hyperbolic tangent]     *
  50.  *                                 *
  51.  *                                 *
  52.  * Description of q-type Data                     *
  53.  * --------------------------                     *
  54.  *                                 *
  55.  * These routines operate on floating point data that consists     *
  56.  * of 24 words of type short per datum.  This data is referred     *
  57.  * to q-type data.  The layout of this data type is:         *
  58.  *                                 *
  59.  *     s   e     o    21-word significand             *
  60.  *     | 1 | 2 | 3 | _____________________ |             *
  61.  *                                 *
  62.  *     s - sign (0 for plus, -1 for minus)             *
  63.  *     e - exponent (0x7fff at most, 0x4001 for 1.0)         *
  64.  *     o - overflow from significand                 *
  65.  *                                 *
  66.  * All routines use pointers to the arrays as arguments.     *
  67.  *                                 *
  68.  * The result is always normalized after each arithmetic     *
  69.  * operation.                             *
  70.  * All arithmetic results are chopped. No rounding is performed  *
  71.  * except on conversion to double precision.             *
  72.  *                                 *
  73.  * ------------------------------------------------------------- */
  74. /* ------------------------------------------------------------- *
  75.  * Revision history:                         *
  76.  *                                 *
  77.  * SLM,  5 Jan 84    PDP-11 assembly language version     *
  78.  * SLM,  2 Mar 86    fixed bug in asctoq()             *
  79.  * SLM,  6 Dec 86    C language version             *
  80.  *                                 *
  81.  * ------------------------------------------------------------- */
  82. #include <stdio.h>
  83. #include "mconf.h"
  84. #include "qhead.h"
  85.  
  86. #ifdef UNK
  87. #if BIGENDIAN
  88. #define MIEEE
  89. #else
  90. #define IBMPC
  91. #endif
  92. #undef UNK
  93. #endif
  94.  
  95. /* array index of the sign (a whole int wasted) */
  96. #define SIGNWORD 0
  97. /* array index of the exponent */
  98. #define E 1
  99. /* array index of first (high order) mantissa word */
  100. #define M 2
  101. /* least significant word */
  102. #define LSWORD NQ-1
  103. /* least significant guard word */
  104. #define LSGUARD NQ
  105. /* Define 1 for sticky bit rounding */
  106. #ifndef STICKY
  107. #define STICKY 0
  108. #endif
  109.  
  110.     /* accumulators */
  111. static QELT ac1[NQ+1] = {0};
  112. static QELT ac2[NQ+1] = {0};
  113. static QELT ac3[NQ+1] = {0};
  114. static QELT ac4[NQ+1] = {0};
  115.     /* shift count register */
  116. static int SC = 0;
  117. extern QELT qzero[], qone[];
  118.  
  119. int shdn1(), shup1(), shdn8(), shup8(), shdn16(), shup16();
  120. int mulm(), mulin(), divm(), qmovz(), normlz();
  121. int cmpm(), addm(), subm(), qadd1(), shift();
  122. int mtherr();
  123.  
  124.  
  125. /* Absolute value  */
  126.  
  127. int qabs (x)
  128. QELT x[];
  129. {
  130.   x[SIGNWORD] = 0;
  131.   return  0;
  132. }
  133.  
  134.  
  135. /* Test if negative  */
  136.  
  137. int qisneg(x)
  138. QELT x[];
  139. {
  140. return (x[SIGNWORD] != 0);
  141. }
  142.  
  143. /*
  144. ;    Negate
  145. */
  146.  
  147. int qneg(x)
  148. QELT x[];
  149. {
  150.  
  151. #if WORDSIZE == 32
  152. x[SIGNWORD] ^= -1; /* complement the sign */
  153. #else
  154. x[SIGNWORD] ^= 0xffff;
  155. #endif
  156. return 0;
  157. }
  158.  
  159. /*
  160. ; convert long integer to q type
  161. ;
  162. ;    long l;
  163. ;    QELT q[NQ];
  164. ;    ltoq( &l, q );
  165. ; note &l is the memory address of l
  166. */
  167.  
  168. int ltoq( lp, y )
  169. long *lp;
  170. QELT y[];
  171. {
  172. long ll;
  173.  
  174. qclear( ac1 );
  175. ll = *lp;        /* local copy of lp */
  176. if( ll < 0 )
  177.     {
  178.     ll = -ll;    /* make it positive */
  179.     ac1[SIGNWORD] = -1; /* put correct sign in the q type number */
  180.     }
  181.  
  182. /* move the long integer to ac1 mantissa area */
  183. #if WORDSIZE == 32
  184. ac1[M] = ll;
  185. ac1[E] = EXPONE - 1;    /* exponent if normalize shift count were 0 */
  186. #else
  187. ac1[M] = ll >> 16;
  188. ac1[M+1] = ll;
  189. ac1[E] = EXPONE+15;
  190. #endif
  191.  
  192. ac1[LSGUARD] = 0;
  193. if( normlz( ac1 ) )    /* normalize the mantissa */
  194.     qclear( ac1 );    /* it was zero */
  195. else
  196.     ac1[E] -= SC;    /* else subtract shift count from exponent */
  197. qmov( ac1, y );        /* output the answer */
  198. return 0;
  199. }
  200.  
  201. /*
  202. ;    convert DEC double precision to q type
  203. ;    double d;
  204. ;    QELT q[NQ];
  205. ;    dtoq( &d, q );
  206. */
  207. int dtoq( d, y )
  208. unsigned short *d;
  209. QELT *y;
  210. {
  211. register QELT r, *p;
  212.  
  213. qclear(y);        /* start with a zero */
  214. p = y;            /* point to our number */
  215. r = *d;            /* get DEC exponent word */
  216. if( r & 0x8000 )
  217.     *p = -1;    /* fill in our sign */
  218. ++p;            /* bump pointer to our exponent word */
  219. r &= 0x7fff;        /* strip the sign bit */
  220. if( r == 0 )        /* answer = 0 if high order DEC word = 0 */
  221.     return 0;
  222.  
  223.  
  224. r >>= 7;    /* shift exponent word down 7 bits */
  225. r -= 0200;    /* subtract DEC exponent offset */
  226. r += EXPONE - 1;    /* add our q type exponent offset */
  227. *p++ = r;    /* to form our exponent */
  228.  
  229. r = *d++;    /* now do the high order mantissa */
  230. r &= 0177;    /* strip off the DEC exponent and sign bits */
  231. r |= 0200;    /* the DEC understood high order mantissa bit */
  232.     /* put result in our high guard word, combined with more mantissa */
  233. #if WORDSIZE == 32
  234. r = (r << 16) | *d++;
  235. *p++ = r;
  236. r = *d++;
  237. r = (r << 16) | *d++;
  238. *p = r;
  239. #else
  240. *p++ = r;    /* put result in our high guard word */
  241. *p++ = *d++;    /* fill in the rest of our mantissa */
  242. *p++ = *d++;
  243. *p = *d;
  244. #endif
  245. shdn8(y);    /* shift our mantissa down 8 bits */
  246. return 0;
  247. }
  248.  
  249. /*
  250. ;    convert q type to DEC double precision
  251. ;    double d;
  252. ;    QELT q[NQ];
  253. ;    qtod( q, &d );
  254. */
  255.  
  256. int qtod( x, d )
  257. QELT *x;
  258. unsigned short *d;
  259. {
  260. register int r;
  261. int i, j;
  262.  
  263. *d = 0;
  264. if( *x != 0 )
  265.     *d = 0100000;
  266. qmovz( x, ac1 );
  267. r = ac1[E];
  268. if( r < (EXPONE - 1 - 0200) )
  269.     goto zout;
  270. #if WORDSIZE == 32
  271. i = ac1[M+2];
  272. #else
  273. i = ac1[M+4];
  274. #endif
  275. if( (i & 0200) != 0 )
  276.     {
  277.     if( (i & 0377) == 0200 )
  278.         {
  279.         if( (i & 0400) != 0 )
  280.             {
  281.         /* check all less significant bits */
  282. #if WORDSIZE == 32
  283.             for( j=M+3; j<=NQ; j++ )
  284. #else
  285.             for( j=M+5; j<=NQ; j++ )
  286. #endif
  287.                 {
  288.                 if( ac1[j] != 0 )
  289.                     goto yesrnd;
  290.                 }
  291.             }
  292.         goto nornd;
  293.         }
  294. yesrnd:
  295.     qclear( ac2 );
  296. #if WORDSIZE == 32
  297.     ac2[ M+2 ] = 0200;
  298. #else
  299.     ac2[ M+4 ] = 0200;
  300. #endif
  301.     ac2[NQ] = 0;
  302.     addm( ac2, ac1 );
  303.     normlz(ac1);
  304.     r -= SC;
  305.     }
  306.  
  307. nornd:
  308.  
  309. r -= EXPONE - 1;
  310. r += 0200;
  311. if( (int) r < 0 )
  312.     {
  313. zout:
  314.     *d++ = 0;
  315.     *d++ = 0;
  316.     *d++ = 0;
  317.     *d++ = 0;
  318.     return 0;
  319.     }
  320. if( r >= 0377 )
  321.     {
  322.     *d++ = 077777;
  323.     *d++ = -1;
  324.     *d++ = -1;
  325.     *d++ = -1;
  326.     return 0;
  327.     }
  328. r &= 0377;
  329. r <<= 7;
  330. shup8( ac1 );
  331. ac1[M] &= 0177;
  332. r |= ac1[M];
  333. *d++ |= r;
  334. #if WORDSIZE == 32
  335. *d++ = ac1[M+1] >> 16;
  336. *d++ = ac1[M+1];
  337. *d++ = ac1[M+2] >> 16;
  338. #else
  339. *d++ = ac1[M+1];
  340. *d++ = ac1[M+2];
  341. *d++ = ac1[M+3];
  342. #endif
  343. return 0;
  344. }
  345.  
  346. /*
  347. ;    Find integer and fractional parts
  348.  
  349. ;    long i;
  350. ;    QELT x[NQ], frac[NQ];
  351. ;    qifrac( x, &i, frac );
  352. */
  353.  
  354. int qifrac( x, i, frac )
  355. QELT x[];
  356. long *i;
  357. QELT frac[];
  358. {
  359.  
  360. qmovz( x, ac1 );
  361. SC = ac1[E] - (EXPONE - 1);
  362. if( SC <= 0 )
  363.     {
  364. /* if exponent <= 0, integer = 0 and argument is fraction */
  365.     *i = 0L;
  366.     qmov( ac1, frac );
  367.     return 0;
  368.     }
  369. if( SC > 31 )
  370.     {
  371. /*
  372. ;    long integer overflow: output large integer
  373. ;    and correct fraction
  374. */
  375.     *i = 0x7fffffff;
  376.     shift( ac1 );
  377.     goto lab10;
  378.     }
  379. #if WORDSIZE == 16
  380. if( SC > 16 )
  381.     {
  382. /*
  383. ; shift more than 16 bits: shift up SC-16, output the integer,
  384. ; then complete the shift to get the fraction.
  385. */
  386.     SC -= 16;
  387.     shift( ac1 );
  388.     *i = ((unsigned long )ac1[M] << 16) | (unsigned short )ac1[M+1];
  389. /*
  390.  *    p = (short *)i;
  391.  * #ifdef DEC
  392.  *    *p++ = ac1[M];
  393.  *    *p++ = ac1[M+1];
  394.  * #else
  395.  *    *p++ = ac1[M+1];
  396.  *    *p++ = ac1[M];
  397.  * #endif
  398.  */
  399.     shup16( ac1 );
  400.     goto lab10;
  401.     }
  402. #endif
  403.  
  404. shift( ac1 );
  405. #if WORDSIZE == 32
  406. *i = ac1[M];
  407. #else
  408. *i = ac1[M] & 0xffff;
  409. #endif
  410.  
  411. lab10:
  412. if( x[SIGNWORD] )
  413.     *i = -(*i);
  414. ac1[SIGNWORD] = 0;
  415. ac1[E] = EXPONE - 1;
  416. ac1[M] = 0;
  417. if( normlz(ac1) )
  418.     qclear( ac1 );
  419. else
  420.     ac1[E] -= SC;
  421.  
  422. qmov( ac1, frac );
  423. return 0;
  424. }
  425.  
  426. int qldexp( x, n )
  427. QELT x[];
  428. int n;
  429. {
  430. int kx, k;
  431.  
  432. kx = x[E];
  433. k = kx    +  n;
  434. if( k < 0 )
  435.   {
  436.     if( n > 0 )
  437.       qinfin(x);
  438.     else
  439.       qclear(x);
  440.     return 0;
  441.   }
  442. x[E] = k;
  443. return 0;
  444. }
  445.  
  446. /*
  447. ;    subtract
  448. ;
  449. ;    QELT a[NQ], b[NQ], c[NQ];
  450. ;    qsub( a, b, c );     c = b - a
  451. */
  452.  
  453. static short subflg = 0;
  454.  
  455. int qsub( a, b, c )
  456. QELT *a, *b, *c;
  457. {
  458.  
  459. subflg = 1;
  460. qadd1( a, b, c );
  461. return 0;
  462. }
  463.  
  464.  
  465. /*
  466. ;    add
  467. ;
  468. ;    QELT a[NQ], b[NQ], c[NQ];
  469. ;    qadd( a, b, c );     c = b + a
  470. */
  471.  
  472. int qadd( a, b, c )
  473. QELT *a, *b, *c;
  474. {
  475.  
  476. subflg = 0;
  477. qadd1( a, b, c );
  478. return 0;
  479. }
  480.  
  481. int qadd1( a, b, c )
  482. QELT *a, *b, *c;
  483. {
  484. long lt;
  485. int i;
  486. #if STICKY
  487. int lost;
  488. #endif
  489.  
  490. qmovz( a, ac1 );
  491. qmovz( b, ac2 );
  492. if( subflg )
  493.     ac1[SIGNWORD] = ~ac1[SIGNWORD];
  494.  
  495. /* compare exponents */
  496. SC = ac1[E] - ac2[E];
  497. if( SC > 0 )
  498.     {    /* put the larger number in ac2 */
  499.     qmovz( ac2, ac3 );
  500.     qmov( ac1, ac2 );
  501.     qmov( ac3, ac1 );
  502.     SC = -SC;
  503.     }
  504.  
  505. #if STICKY
  506. lost = 0;
  507. #endif
  508. if( SC != 0 )
  509.     {
  510.     if( SC < -NBITS-1 )
  511.         goto done;    /* answer same as larger addend */
  512. #if STICKY
  513.     lost = shift( ac1, SC ); /* shift the smaller number down */
  514. #else
  515.     shift( ac1, SC ); /* shift the smaller number down */
  516. #endif
  517.     }
  518. else
  519.     {
  520. /* exponents were the same, so must compare mantissae */
  521.     i = cmpm( ac1, ac2 );
  522.     if( i == 0 )
  523.         {    /* the numbers are identical */
  524.         /* if different signs, result is zero */
  525.         if( ac1[SIGNWORD] != ac2[SIGNWORD] )
  526.             goto underf;
  527.         /* if exponents zero, result is zero */
  528.         if( ac1[E] == 0 )
  529.             goto underf;
  530.         /* if same sign, result is double */
  531.         ac2[E] += 1;
  532.         goto done;
  533.         }
  534.     if( i > 0 )
  535.         {    /* put the larger number in ac2 */
  536.         qmovz( ac2, ac3 );
  537.         qmov( ac1, ac2 );
  538.         qmov( ac3, ac1 );
  539.         }
  540.     }
  541.  
  542. if( ac1[SIGNWORD] == ac2[SIGNWORD] )
  543.     {
  544.     addm( ac1, ac2 );
  545.     subflg = 0;
  546.     }
  547. else
  548.     {
  549.     subm( ac1, ac2 );
  550.     subflg = 1;
  551.     }
  552. if( normlz(ac2) )
  553.     goto underf;
  554.  
  555. lt = (long )ac2[E] - SC;
  556. if( lt > MAXEXP )
  557.     goto overf;
  558. if( lt < 0 )
  559.     {
  560. /*    mtherr( "qadd", UNDERFLOW );*/
  561.     goto underf;
  562.     }
  563. ac2[E] = lt;
  564.  
  565. /* round off */
  566. #if WORDSIZE == 32
  567. i = ac2[LSGUARD];
  568. #else
  569. i = ac2[NQ] & 0xffff;
  570. #endif
  571. if( i & SIGNBIT )
  572.     {
  573. #if STICKY
  574.     if( i == SIGNBIT )
  575.         {
  576.         if( lost == 0 )
  577.             {
  578. /* Critical case, round to even */
  579.             if( (ac2[LSWORD] & 1) == 0 )
  580.                 goto done;
  581.             }
  582.         else
  583.             {
  584.             if( subflg != 0 )
  585.                 goto done;
  586.             }
  587.         }
  588. #else
  589.     if( subflg != 0 )
  590.         goto done;
  591. #endif
  592.  
  593.     qclear( ac1 );
  594.     ac1[LSGUARD] = 0;
  595.     ac1[LSWORD] = 1;
  596.     addm( ac1, ac2 );
  597.     normlz(ac2);
  598.     if( SC )
  599.         {
  600.         lt = (long )ac2[E] - SC;
  601.         if( lt > MAXEXP )
  602.         goto overf;
  603.         ac2[E] = lt;
  604.         }
  605.     }
  606. done:
  607. qmov( ac2, c );
  608. return 0;
  609.  
  610. underf:
  611. qclear(c);
  612. return 0;
  613.  
  614. overf:
  615. mtherr( "qadd", OVERFLOW );
  616. qinfin(c);
  617. return 0;
  618. }
  619.  
  620. /*
  621. ;    divide
  622. ;
  623. ;    QELT a[NQ], b[NQ], c[NQ];
  624. ;    qdiv( a, b, c );    c = b / a
  625. */
  626. /* for Newton iteration version:
  627.  * extern short qtwo[];
  628.  * static short qt[NQ] = {0};
  629.  * static short qu[NQ] = {0};
  630.  */
  631. int qdiv( a, b, c )
  632. QELT *a, *b, *c;
  633. {
  634. long lt;
  635.  
  636. if( b[E] == 0 )
  637.     {
  638.     qclear(c);    /* numerator is zero */
  639.     return 0;
  640.     }
  641.  
  642. if( a[E] == 0 )
  643.     {    /* divide by zero */
  644.     qinfin(c);
  645.     mtherr( "qdiv", SING );
  646.     return 0;
  647.     }
  648. qmovz( b, ac3 );
  649. divm( a, ac3 );
  650. /* calculate exponent */
  651. lt = (long )ac3[E] - (long )a[E];
  652. ac3[E] = lt;
  653. ac3[LSGUARD] = 0;
  654. normlz(ac3);
  655. /* Newton iteration version */
  656. lt = lt - SC + EXPONE + 1;
  657. /* Taylor series version */
  658. /*lt = lt - SC + 16385L;*/
  659. ac3[E] = lt;
  660.  
  661. if( a[SIGNWORD] == b[SIGNWORD] )
  662.     ac3[SIGNWORD] = 0;
  663. else
  664.     ac3[SIGNWORD] = -1;
  665.  
  666. qmov( ac3, c );
  667. return 0;
  668. }
  669.  
  670. /*
  671. ;    multiply
  672. ;
  673. ;    QELT a[NQ], b[NQ], c[NQ];
  674. ;    qmul( a, b, c );    c = b * a
  675. */
  676. int qmul( a, b, c )
  677. QELT *a, *b, *c;
  678. {
  679. QELT *p;
  680. register int ctr;
  681. long lt;
  682.  
  683. if( (a[E] == 0) || (b[E] == 0) )
  684.     {
  685.     qclear(c);
  686.     return 0;
  687.     }
  688. /* detect multiplication by small integer a */
  689. if( a[M+2] == 0 )
  690.     {
  691.     p = &a[M+3];
  692.     for( ctr=M+3; ctr<NQ; ctr++ )
  693.         {
  694.         if( *p++ != 0 )
  695.             goto nota;
  696.         }
  697.     qmov( b, ac3 );
  698.     mulin( a, ac3 );
  699.     lt = ((long)a[E] - (EXPONE-1)) + ((long )ac3[E] - (EXPONE - 1));
  700.     goto mulcon;
  701.     }
  702.  
  703. nota:
  704. /* detect multiplication by small integer b */
  705. if( b[M+2] == 0 )
  706.     {
  707.     p = &b[M+3];
  708.     for( ctr=M+3; ctr<NQ; ctr++ )
  709.         {
  710.         if( *p++ != 0 )
  711.             goto notb;
  712.         }
  713.     qmov( a, ac3 );
  714.     mulin( b, ac3 );
  715.     lt = ((long)b[E] - (EXPONE-1)) + ((long )ac3[E] - (EXPONE - 1));
  716.     goto mulcon;
  717.     }
  718.  
  719. notb:
  720.  
  721. qmov( a, ac3 );
  722. mulm( b, ac3 );
  723. lt = ((long)b[E] - (EXPONE-1)) + ((long )ac3[E] - (EXPONE - 1));
  724.  
  725. mulcon:
  726. /* calculate sign of product */
  727. if( b[SIGNWORD] == a[SIGNWORD] )
  728.     ac3[SIGNWORD] = 0;
  729. else
  730.     ac3[SIGNWORD] = -1;
  731.  
  732. if( lt > MAXEXP )
  733.     goto overf;
  734. ac3[E] = lt;
  735. ac3[NQ] = 0;
  736. if( normlz(ac3) )
  737.     goto underf;
  738. lt = lt - SC + EXPONE -1;
  739. if( lt > MAXEXP )
  740.     goto overf;
  741. if( lt < 0 )
  742.     goto underf;
  743. ac3[E] = lt;
  744. qmov( ac3, c );
  745. return 0;
  746.  
  747. underf:
  748. qclear(c);
  749. return 0;
  750.  
  751. overf:
  752. qinfin(c);
  753. mtherr( "qmul", OVERFLOW );
  754. return 0;
  755. }
  756.  
  757.  
  758.  
  759.  
  760. /* Multiply, a has at most WORDSIZE significant bits */
  761.  
  762. int qmuli( a, b, c )
  763. QELT *a, *b, *c;
  764. {
  765. long lt;
  766.  
  767. if( (a[E] == 0) || (b[E] == 0) )
  768.     {
  769.     qclear(c);
  770.     return 0;
  771.     }
  772.  
  773. qmov( b, ac3 );
  774. mulin( a, ac3 );
  775.  
  776. /* calculate sign of product */
  777. if( b[SIGNWORD] == a[SIGNWORD] )
  778.     ac3[SIGNWORD] = 0;
  779. else
  780.     ac3[SIGNWORD] = -1;
  781.  
  782. /* calculate exponent */
  783. lt = ((long)ac3[E] - (EXPONE-1)) + ((long )a[E] - (EXPONE - 1));
  784. if( lt > MAXEXP )
  785.     goto overf;
  786. ac3[E] = lt;
  787. ac3[NQ] = 0;
  788. if( normlz(ac3) )
  789.     goto underf;
  790. lt = lt - SC + EXPONE - 1;
  791. if( lt > MAXEXP )
  792.     goto overf;
  793. if( lt < 0 )
  794.     goto underf;
  795. ac3[E] = lt;
  796. qmov( ac3, c );
  797. return 0;
  798.  
  799. underf:
  800. qclear(c);
  801. return 0;
  802.  
  803. overf:
  804. qinfin(c);
  805. mtherr( "qmuli", OVERFLOW );
  806. return 0;
  807. }
  808.  
  809.  
  810.  
  811.  
  812. /*
  813. ;    Compare mantissas
  814. ;
  815. ;    QELT a[NQ], b[NQ];
  816. ;    cmpm( a, b );
  817. ;
  818. ;    for the mantissas:
  819. ;    returns    +1 if a > b
  820. ;         0 if a == b
  821. ;        -1 if a < b
  822. */
  823.  
  824. int cmpm( a, b )
  825. register QELT *a, *b;
  826. {
  827. int i;
  828.  
  829. a += M; /* skip up to mantissa area */
  830. b += M;
  831. for( i=0; i<OMG; i++ )
  832.     {
  833.     if( *a++ != *b++ )
  834.         goto difrnt;
  835.     }
  836. return(0);
  837.  
  838. difrnt:
  839. if( (unsigned int) *(--a) > (unsigned int) *(--b) )
  840.     return(1);
  841. else
  842.     return(-1);
  843. }
  844.  
  845. /*
  846. ;    shift mantissa
  847. ;
  848. ;    Shifts mantissa area up or down by the number of bits
  849. ;    given by the variable SC.
  850. */
  851.  
  852. int shift( x )
  853. QELT *x;
  854. {
  855. QELT *p;
  856. #if STICKY
  857. int lost;
  858. #endif
  859.  
  860. if( SC == 0 )
  861.     return(0);
  862.  
  863. #if STICKY
  864. lost = 0;
  865. #endif
  866. if( SC < 0 )
  867.     {
  868.     p = x + LSGUARD;
  869.     SC = -SC;
  870.     while( SC >= 16 )
  871.         {
  872. #if STICKY
  873.         lost |= *p;
  874. #endif
  875.         shdn16(x);
  876.         SC -= 16;
  877.         }
  878.  
  879.     while( SC >= 8 )
  880.         {
  881. #if STICKY
  882.         lost |= *p & 0xff;
  883. #endif
  884.         shdn8(x);
  885.         SC -= 8;
  886.         }
  887.  
  888.     while( SC > 0 )
  889.         {
  890. #if STICKY
  891.         lost |= *p & 1;
  892. #endif
  893.         shdn1(x);
  894.         SC -= 1;
  895.         }
  896.     }
  897. else
  898.     {
  899.     while( SC >= 16 )
  900.         {
  901.         shup16(x);
  902.         SC -= 16;
  903.         }
  904.  
  905.     while( SC >= 8 )
  906.         {
  907.         shup8(x);
  908.         SC -= 8;
  909.         }
  910.  
  911.     while( SC > 0 )
  912.         {
  913.         shup1(x);
  914.         SC -= 1;
  915.         }
  916.     }
  917. #if STICKY
  918. return( lost );
  919. #else
  920. return(0);
  921. #endif
  922. }
  923.  
  924. /*
  925. ;    normalize
  926. ;
  927. ; shift normalizes the mantissa area pointed to by R1
  928. ; shift count (up = positive) returned in SC
  929. */
  930.  
  931. int normlz(x)
  932. QELT x[];
  933. {
  934. register QELT *p;
  935.  
  936. SC = 0;
  937. p = &x[M];
  938. if( *p != 0 )
  939.     goto normdn;
  940. ++p;
  941. if( *p & SIGNBIT )
  942.     return(0);    /* already normalized */
  943.  
  944. while( *p == 0 )
  945.     {
  946.     shup16(x);
  947.     SC += 16;
  948. /* With guard word, there are NBITS+WORDSIZE bits available.
  949.  * return true if all are zero.
  950.  */
  951.     if( SC > NBITS )
  952.         return(1);
  953.     }
  954.  
  955. /* see if high byte is zero */
  956. #if WORDSIZE == 32
  957. while( (*p & 0xff000000) == 0 )
  958. #else
  959. while( (*p & 0xff00) == 0 )
  960. #endif
  961.     {
  962.     shup8(x);
  963.     SC += 8;
  964.     }
  965.  
  966. /* now shift 1 bit at a time */
  967. while( (*p & SIGNBIT) == 0)
  968.     {
  969.     shup1(x);
  970.     SC += 1;
  971. /*
  972.     if( SC > NBITS )
  973.         {
  974.         printf( "normlz error\n");
  975.         return(0);
  976.         }
  977. */
  978.     }
  979. return(0);
  980.  
  981. /* normalize by shifting down out of the high guard word
  982.    of the mantissa */
  983.  
  984. normdn:
  985.  
  986. #if WORDSIZE == 32
  987. while( *p & 0xffffff00 )
  988. #else
  989. while( *p & 0xff00 )
  990. #endif
  991.     {
  992.     shdn8(x);
  993.     SC -= 8;
  994.     }
  995. while( *p != 0 )
  996.     {
  997.     shdn1(x);
  998.     SC -= 1;
  999. /*
  1000.     if( SC < -NBITS )
  1001.         {
  1002.         printf( "low normlz error\n");
  1003.         return(0);
  1004.         }
  1005. */
  1006.     }
  1007. return(0);
  1008. }
  1009.  
  1010. /*
  1011. ; Clear out entire number, including sign and exponent, pointed
  1012. ; to by x
  1013. ;
  1014. ; QELT x[];
  1015. ; qclear( x );
  1016. */
  1017. /* Moved to qfltb.c */
  1018. #if 0
  1019. int qclear( x )
  1020. register QELT *x;
  1021. {
  1022. register int i;
  1023.  
  1024. for( i=0; i<NQ; i++ )
  1025.     *x++ = 0;
  1026. return 0;
  1027. }
  1028. #endif
  1029.  
  1030. /*
  1031. ; Fill entire number, including exponent and mantissa, with
  1032. ; largest possible number.
  1033. */
  1034.  
  1035. int qinfin(x)
  1036. QELT *x;
  1037. {
  1038. register int i;
  1039.  
  1040. ++x;    /* skip over the sign */
  1041. *x++ = MAXEXP;
  1042. *x++ = 0;
  1043. for( i=0; i<OMG-1; i++ )
  1044.     *x++ = -1;
  1045. return 0;
  1046. }
  1047.  
  1048. /* normalization program */
  1049. int qnrmlz(x)
  1050. QELT *x;
  1051. {
  1052.  
  1053. qmovz( x, ac1 );
  1054. normlz( ac1 );    /* shift normalize the mantissa */
  1055. ac1[E] -= SC;    /* subtract the shift count from the exponent */
  1056. qmov( ac1, x );
  1057. return 0;
  1058. }
  1059.  
  1060.  
  1061.  
  1062.  
  1063. /*
  1064. ; Convert IEEE double precision to Q type
  1065. ;    double d;
  1066. ;    QELT q[N+2];
  1067. ;    etoq( &d, q );
  1068. */
  1069.  
  1070. int etoq( e, y )
  1071. unsigned short *e;
  1072. QELT *y;
  1073. {
  1074. #ifdef DEC
  1075. dtoq(e,y);
  1076. #else
  1077. register int r;
  1078. register QELT *p;
  1079. QELT yy[NQ+1];
  1080. int denorm;
  1081.  
  1082.  
  1083. denorm = 0;    /* flag if denormalized number */
  1084. qclear(yy);
  1085. yy[NQ] = 0;
  1086.  
  1087. #ifdef MIEEE
  1088.  
  1089. #endif
  1090.  
  1091. #ifdef IBMPC
  1092. e += M+1;
  1093. #endif
  1094.  
  1095. /*
  1096. r = *e & 0x7fff;
  1097. if( r == 0 )
  1098.     return 0;
  1099. */
  1100.  
  1101. r = *e;
  1102. yy[SIGNWORD] = 0;
  1103. if( r & 0x8000 )
  1104.     yy[SIGNWORD] = -1;
  1105.  
  1106. yy[M] = (r & 0x0f) | 0x10;
  1107. r &= ~0x800f;    /* strip sign and 4 mantissa bits */
  1108. r >>= 4;
  1109. /* If zero exponent, then the mantissa is denormalized.
  1110.  * So take back the understood high mantissa bit. */
  1111. if( r == 0 )
  1112.     {
  1113.     denorm = 1;
  1114.     yy[M] &= ~0x10;
  1115.     }
  1116. r += EXPONE - 01777;
  1117. yy[E] = r;
  1118. p = &yy[M+1];
  1119. #if WORDSIZE == 32
  1120. #ifdef MIEEE
  1121. ++e;
  1122. *p = ((unsigned int) *e++) << 16;
  1123. *p++ |= *e++;
  1124. *p = ((unsigned int) *e++) << 16;
  1125. #endif
  1126. #ifdef IBMPC
  1127. *p = ((unsigned int) *(--e)) << 16;
  1128. *p++ |= *(--e);
  1129. *p = ((unsigned int) *(--e)) << 16;
  1130. #endif
  1131. #else
  1132. #ifdef MIEEE
  1133. ++e;
  1134. *p++ = *e++;
  1135. *p++ = *e++;
  1136. *p++ = *e++;
  1137. #endif
  1138. #ifdef IBMPC
  1139. *p++ = *(--e);
  1140. *p++ = *(--e);
  1141. *p++ = *(--e);
  1142. #endif
  1143. #endif
  1144. SC = -5;
  1145. shift(yy);
  1146. if( denorm )
  1147.     { /* if zero exponent, then normalize the mantissa */
  1148.     if( normlz( yy ) )
  1149.         qclear(yy);
  1150.     else
  1151.         yy[E] -= SC-1;
  1152.     }
  1153. qmov( yy, y );
  1154. /* not DEC */
  1155. #endif
  1156. return 0;
  1157. }
  1158.  
  1159. /*
  1160. ; Q type to IEEE double precision
  1161. ;    double d;
  1162. ;    QELT q[N+2];
  1163. ;    qtoe( q, &d );
  1164. */
  1165. int qtoe( x, e )
  1166. QELT *x;
  1167. short *e;
  1168. {
  1169. #ifdef DEC
  1170. qtod(x,e);
  1171. #else
  1172. int i, j, k;
  1173. register QELT *p;
  1174.  
  1175.  
  1176. #ifdef MIEEE
  1177.  
  1178. #endif
  1179. #ifdef IBMPC
  1180. e += M+1;
  1181. #endif
  1182.  
  1183.  
  1184. *e = 0;    /* output high order */
  1185. p = &ac1[SIGNWORD];
  1186. qmovz( x, ac1 );
  1187. if( *p++ != 0 )
  1188.     *e = 0x8000;    /* output sign bit */
  1189.  
  1190. normlz(ac1);
  1191. *p -= SC;
  1192.  
  1193. i = *p++ -(EXPONE - 1023);    /* adjust exponent for offsets */
  1194.  
  1195. /* Handle denormalized small numbers.
  1196.  * The tiniest number represented appears to be 2**-1074.
  1197.  */
  1198. if( i <= 0 )
  1199.     {
  1200.     if( i > -53 )
  1201.         {
  1202.         SC = i - 1;
  1203.         shift( ac1 );
  1204.         i = 0;
  1205.         }
  1206.     else
  1207.         {
  1208. /*ozero:*/
  1209. #ifdef MIEEE
  1210.         ++e;
  1211.         *e++ = 0;
  1212.         *e++ = 0;
  1213.         *e++ = 0;
  1214. #endif
  1215. #ifdef IBMPC
  1216.         *(--e) = 0;
  1217.         *(--e) = 0;
  1218.         *(--e) = 0;
  1219. #endif
  1220.         return 0;
  1221.         }
  1222.     }
  1223.  
  1224. /* round off to nearest or even */
  1225. #if WORDSIZE == 32
  1226. k = ac1[M+2];
  1227. #else
  1228. k = ac1[M+4];
  1229. #endif
  1230. if( (k & 0x400) != 0 )
  1231.     {
  1232.     if( (k & 0x07ff) == 0x400 )
  1233.         {
  1234.         if( (k & 0x800) != 0 )
  1235.             {
  1236.         /* check all less significant bits */
  1237. #if WORDSIZE == 32
  1238.             for( j=M+3; j<=NQ; j++ )
  1239. #else
  1240.             for( j=M+5; j<=NQ; j++ )
  1241. #endif
  1242.                 {
  1243.                 if( ac1[j] != 0 )
  1244.                     goto yesrnd;
  1245.                 }
  1246.             }
  1247.         goto nornd;
  1248.         }
  1249. yesrnd:
  1250.     qclear( ac2 );
  1251.     ac2[NQ] = 0;
  1252. #if WORDSIZE == 32
  1253.     ac2[M+2] = 0x800;
  1254. #else
  1255.     ac2[M+4] = 0x800;
  1256. #endif
  1257.     addm( ac2, ac1 );
  1258.     if( ac1[2] )
  1259.         {
  1260.         shdn1(ac1);
  1261.         i += 1;
  1262.         }
  1263.     if( (i == 0) && (ac1[M+1] & SIGNBIT) )
  1264.         i += 1;
  1265.     }
  1266.  
  1267. nornd:
  1268.  
  1269. if( i > 2047 )
  1270.     {    /* Saturate at largest number less than infinity. */
  1271.     mtherr( "qtoe", OVERFLOW );
  1272.     *e |= 0x7fef;
  1273. #ifdef MIEEE
  1274.     ++e;
  1275.     *e++ = 0xffff;
  1276.     *e++ = 0xffff;
  1277.     *e++ = 0xffff;
  1278. #endif
  1279. #ifdef IBMPC
  1280.     *(--e) = 0xffff;
  1281.     *(--e) = 0xffff;
  1282.     *(--e) = 0xffff;
  1283. #endif
  1284.     return 0;
  1285.     }
  1286.  
  1287.  
  1288. i <<= 4;
  1289. SC = 5;
  1290. shift( ac1 );
  1291. i |= *p++ & 0x0f;    /* *p = ac1[M] */
  1292. *e |= i;    /* high order output already has sign bit set */
  1293. #if WORDSIZE == 32
  1294. #ifdef MIEEE
  1295. ++e;
  1296. *e++ = *p >> 16;
  1297. *e++ = *p++;
  1298. *e++ = *p >> 16;
  1299. #endif
  1300. #ifdef IBMPC
  1301. *(--e) = *p >> 16;
  1302. *(--e) = *p++;
  1303. *(--e) = *p >> 16;
  1304. #endif
  1305. #else
  1306. /* 16 bit words */
  1307. #ifdef MIEEE
  1308. ++e;
  1309. *e++ = *p++;
  1310. *e++ = *p++;
  1311. *e++ = *p++;
  1312. #endif
  1313. #ifdef IBMPC
  1314. *(--e) = *p++;
  1315. *(--e) = *p++;
  1316. *(--e) = *p;
  1317. #endif
  1318. #endif
  1319. /* not DEC */
  1320. #endif
  1321. return 0;
  1322. }
  1323.  
  1324. /*
  1325. ; Convert 80-bit IEEE long double precision to Q type
  1326. ;    long double d;
  1327. ;    QELT q[N+2];
  1328. ;    etoq( &d, q );
  1329. */
  1330.  
  1331. int e64toq( e, y )
  1332. unsigned short *e;
  1333. QELT *y;
  1334. {
  1335. register int r;
  1336. register QELT *p;
  1337. QELT yy[NQ+1];
  1338. int denorm;
  1339.  
  1340.  
  1341. denorm = 0;    /* flag if denormalized number */
  1342. qclear(yy);
  1343. yy[NQ] = 0;
  1344.  
  1345. #ifdef MIEEE
  1346.  
  1347. #endif
  1348.  
  1349.  
  1350. #ifdef IBMPC
  1351. e += 4;
  1352. #endif
  1353.  
  1354. r = *e;
  1355. yy[SIGNWORD] = 0;
  1356. if( r & 0x8000 )
  1357.     yy[SIGNWORD] = -1;
  1358.  
  1359. r &= 0x7fff;
  1360. /* If zero exponent, then the mantissa is denormalized. */
  1361. if( r == 0 )
  1362.     {
  1363.     denorm = 1;
  1364.     }
  1365. r += (EXPONE-0x3fff);
  1366. yy[E] = r;
  1367. p = &yy[M+1];
  1368. #if WORDSIZE== 32
  1369. #ifdef MIEEE
  1370. ++e;
  1371. *p = ((unsigned int) *e++) << 16;
  1372. *p++ |= *e++;
  1373. *p = ((unsigned int) *e++) << 16;
  1374. *p++ |= *e++;
  1375. #endif
  1376. #ifdef IBMPC
  1377. *p = ((unsigned int) *(--e)) << 16;
  1378. *p++ |= *(--e);
  1379. *p = ((unsigned int) *(--e)) << 16;
  1380. *p++ |= *(--e);
  1381. #endif
  1382. #else
  1383. /* 16-bit wordsize */
  1384. #ifdef MIEEE
  1385. ++e;
  1386. *p++ = *e++;
  1387. *p++ = *e++;
  1388. *p++ = *e++;
  1389. *p++ = *e++;
  1390. #endif
  1391. #ifdef IBMPC
  1392. *p++ = *(--e);
  1393. *p++ = *(--e);
  1394. *p++ = *(--e);
  1395. *p++ = *(--e);
  1396. #endif
  1397. #endif
  1398. if( denorm )
  1399.     { /* if zero exponent, then normalize the mantissa */
  1400. #ifdef IBMPC
  1401. /* For Intel long double, shift denormal significand up 1
  1402.    -- but only if the top significand bit is zero.  */
  1403.     if( (yy[M+1] & SIGNBIT) == 0 )
  1404.         shup1( yy );
  1405. #endif
  1406.     if( normlz( yy ) )
  1407.         qclear(yy);
  1408.     else
  1409.         yy[E] -= SC-1;
  1410.     }
  1411. qmov( yy, y );
  1412. return 0;
  1413. }
  1414.  
  1415. /*
  1416. ; Q type to 80-bit IEEE long double precision
  1417. ;    long double d;
  1418. ;    QELT q[N+2];
  1419. ;    qtoe( q, &d );
  1420. */
  1421. int qtoe64( x, e )
  1422. QELT *x;
  1423. unsigned short *e;
  1424. {
  1425. QELT k;
  1426. short i, j;
  1427. register QELT *p;
  1428.  
  1429.  
  1430. #ifdef MIEEE
  1431.  
  1432. #endif
  1433. #ifdef IBMPC
  1434. e += 4;
  1435. #if 0
  1436. /* NOTE: if data type is 96 bits wide, clear the last word here. */
  1437. *(e+1)= 0;
  1438. #endif
  1439. #endif
  1440.  
  1441.  
  1442. *e = 0;    /* output high order */
  1443. p = &ac1[0];
  1444. qmovz( x, p );
  1445. if( ac1[SIGNWORD] )
  1446.     *e = 0x8000;    /* output sign bit */
  1447. ++p;
  1448. normlz(ac1);
  1449. *p -= SC;
  1450.  
  1451. i = *p++ - (EXPONE-0x3fff);    /* adjust exponent for offsets */
  1452.  
  1453. /* We can't handle denormal numbers. */
  1454. if( i <= 0 )
  1455.     {
  1456. /*    if( i > -2 ) */
  1457.     if( i > -64 )
  1458.         {
  1459.         SC = i - 1;
  1460.         shift( ac1 );
  1461.         i = 0;
  1462.         }
  1463.     else
  1464.         {
  1465. /*ozero:*/
  1466. #ifdef MIEEE
  1467.         ++e;
  1468.         *e++ = 0;
  1469.         *e++ = 0;
  1470.         *e++ = 0;
  1471.         *e++ = 0;
  1472. #endif
  1473. #ifdef IBMPC
  1474.         *(--e) = 0;
  1475.         *(--e) = 0;
  1476.         *(--e) = 0;
  1477.         *(--e) = 0;
  1478. #endif
  1479.         return 0;
  1480.         }
  1481.     }
  1482.  
  1483. /* round off to nearest or even */
  1484. #if WORDSIZE == 32
  1485. k = ac1[M+3];
  1486. #else
  1487. k = ac1[M+5];
  1488. #endif
  1489. if( (k & SIGNBIT) != 0 )
  1490.     {
  1491.     if( (k & ((QELT) -1)) == SIGNBIT )
  1492.         {
  1493.         /* check all less significant bits */
  1494. #if WORDSIZE == 32
  1495. /* 0 1 2  3   4   5   6 */
  1496. /* S E M 1,2 3,4 5,6 7,8 */
  1497.         for( j=M+4; j<=NQ; j++ )
  1498. #else
  1499. /* 0 1 2 3 4 5 6 7 8 9 10 */
  1500. /* S E M 1 2 3 4 5 6 7 8 */
  1501.         for( j=M+6; j<=NQ; j++ )
  1502. #endif
  1503.             {
  1504.             if( ac1[j] != 0 )
  1505.                 goto yesrnd;
  1506.             }
  1507. /* round to even */
  1508. #if WORDSIZE == 32
  1509.         if( (ac1[4] & 1) == 0 )
  1510. #else
  1511.         if( (ac1[6] & 1) == 0 )
  1512. #endif
  1513.           goto nornd;
  1514.         }
  1515. yesrnd:
  1516.     qclear( ac2 );
  1517.     ac2[NQ] = 0;
  1518. #if WORDSIZE == 32
  1519.     ac2[M+3] = SIGNBIT;
  1520. #else
  1521.     ac2[M+5] = SIGNBIT;
  1522. #endif
  1523.     addm( ac2, ac1 );
  1524.     if( ac1[M] )
  1525.         {
  1526.         shdn1(ac1);
  1527.         i += 1;
  1528.         }
  1529.     if( (i == 0) && (ac1[3] & SIGNBIT) )
  1530.         i += 1;
  1531.     }
  1532.  
  1533. nornd:
  1534.  
  1535. *e |= i;    /* high order output already has sign bit set */
  1536. p = &ac1[M+1];
  1537. #if WORDSIZE == 32
  1538. #ifdef MIEEE
  1539. ++e;
  1540. *e++ = *p >> 16;
  1541. *e++ = *p++;
  1542. *e++ = *p >> 16;
  1543. *e++ = *p++;
  1544. #endif
  1545. #ifdef IBMPC
  1546. *(--e) = *p >> 16;
  1547. *(--e) = *p++;
  1548. *(--e) = *p >> 16;
  1549. *(--e) = *p++;
  1550. #endif
  1551. #else
  1552. /* 16-bit words */
  1553. #ifdef MIEEE
  1554. ++e;
  1555. *e++ = *p++;
  1556. *e++ = *p++;
  1557. *e++ = *p++;
  1558. *e++ = *p;
  1559. #endif
  1560. #ifdef IBMPC
  1561. *(--e) = *p++;
  1562. *(--e) = *p++;
  1563. *(--e) = *p++;
  1564. *(--e) = *p;
  1565. #endif
  1566. #endif
  1567. return 0;
  1568. }
  1569.  
  1570. /*
  1571. ; Convert 128-bit IEEE long double precision to Q type
  1572. ;    long double d;
  1573. ;    QELT q[N+2];
  1574. ;    etoq( &d, q );
  1575. */
  1576.  
  1577. int e113toq( e, y )
  1578. unsigned short *e;
  1579. QELT *y;
  1580. {
  1581. register int r;
  1582. register QELT *p;
  1583. QELT yy[NQ+1];
  1584. int denorm;
  1585.  
  1586.  
  1587. denorm = 0;    /* flag if denormalized number */
  1588. qclear(yy);
  1589. yy[NQ] = 0;
  1590.  
  1591. #ifdef MIEEE
  1592.  
  1593. #endif
  1594.  
  1595.  
  1596. #ifdef IBMPC
  1597. e += 7;
  1598. #endif
  1599.  
  1600. r = *e;
  1601. yy[SIGNWORD] = 0;
  1602. if( r & 0x8000 )
  1603.     yy[SIGNWORD] = -1;
  1604.  
  1605. r &= 0x7fff;
  1606. yy[M] = 1;
  1607. /* If zero exponent, then the mantissa is denormalized.
  1608.  * So take back the understood high mantissa bit. */
  1609. if( r == 0 )
  1610.     {
  1611.     denorm = 1;
  1612.     yy[M] = 0;
  1613.     }
  1614. r += 2;
  1615. yy[E] = r;
  1616. p = &yy[M+1];
  1617. #if WORDSIZE== 32
  1618. #ifdef MIEEE
  1619. ++e;
  1620. *p = ((unsigned int) *e++) << 16;
  1621. *p++ |= *e++;
  1622. *p = ((unsigned int) *e++) << 16;
  1623. *p++ |= *e++;
  1624. *p = ((unsigned int) *e++) << 16;
  1625. *p++ |= *e++;
  1626. *p = ((unsigned int) *e++) << 16;
  1627. #endif
  1628. #ifdef IBMPC
  1629. *p = ((unsigned int) *(--e)) << 16;
  1630. *p++ |= *(--e);
  1631. *p = ((unsigned int) *(--e)) << 16;
  1632. *p++ |= *(--e);
  1633. *p = ((unsigned int) *(--e)) << 16;
  1634. *p++ |= *(--e);
  1635. *p = ((unsigned int) *(--e)) << 16;
  1636. #endif
  1637. #else
  1638. /* 16-bit wordsize */
  1639. #ifdef MIEEE
  1640. ++e;
  1641. *p++ = *e++;
  1642. *p++ = *e++;
  1643. *p++ = *e++;
  1644. *p++ = *e++;
  1645. *p++ = *e++;
  1646. *p++ = *e++;
  1647. *p++ = *e++;
  1648. #endif
  1649. #ifdef IBMPC
  1650. *p++ = *(--e);
  1651. *p++ = *(--e);
  1652. *p++ = *(--e);
  1653. *p++ = *(--e);
  1654. *p++ = *(--e);
  1655. *p++ = *(--e);
  1656. *p++ = *(--e);
  1657. #endif
  1658. #endif
  1659. SC = -1;
  1660. shift(yy);
  1661. if( denorm )
  1662.     { /* if zero exponent, then normalize the mantissa */
  1663.     if( normlz( yy ) )
  1664.         qclear(yy);
  1665.     else
  1666.         yy[E] -= SC-1;
  1667.     }
  1668. qmov( yy, y );
  1669. return 0;
  1670. }
  1671.  
  1672. /*
  1673. ; Q type to 128-bit IEEE long double precision
  1674. ;    long double d;
  1675. ;    QELT q[N+2];
  1676. ;    qtoe( q, &d );
  1677. */
  1678. int qtoe113( x, e )
  1679. QELT *x;
  1680. unsigned short *e;
  1681. {
  1682. QELT k;
  1683. short i, j;
  1684. register QELT *p;
  1685.  
  1686.  
  1687. #ifdef MIEEE
  1688.  
  1689. #endif
  1690. #ifdef IBMPC
  1691. e += 7;
  1692. #endif
  1693.  
  1694.  
  1695. *e = 0;    /* output high order */
  1696. p = &ac1[0];
  1697. qmovz( x, p );
  1698. if( ac1[SIGNWORD] )
  1699.     *e = 0x8000;    /* output sign bit */
  1700. ++p;
  1701. normlz(ac1);
  1702. *p -= SC;
  1703.  
  1704. i = *p++ - 2;    /* adjust exponent for offsets */
  1705.  
  1706. /* We can't handle denormal numbers. */
  1707. if( i <= 0 )
  1708.     {
  1709.     if( i > -2 )
  1710.         {
  1711.         SC = i - 1;
  1712.         shift( ac1 );
  1713.         i = 0;
  1714.         }
  1715.     else
  1716.         {
  1717. /*ozero:*/
  1718. #ifdef MIEEE
  1719.         ++e;
  1720.         *e++ = 0;
  1721.         *e++ = 0;
  1722.         *e++ = 0;
  1723.         *e++ = 0;
  1724.         *e++ = 0;
  1725.         *e++ = 0;
  1726.         *e++ = 0;
  1727. #endif
  1728. #ifdef IBMPC
  1729.         *(--e) = 0;
  1730.         *(--e) = 0;
  1731.         *(--e) = 0;
  1732.         *(--e) = 0;
  1733.         *(--e) = 0;
  1734.         *(--e) = 0;
  1735.         *(--e) = 0;
  1736. #endif
  1737.         return 0;
  1738.         }
  1739.     }
  1740.  
  1741. /* round off to nearest or even */
  1742. #if WORDSIZE == 32
  1743. k = ac1[M+4];
  1744. #else
  1745. k = ac1[M+8];
  1746. #endif
  1747. if( (k & 0x4000) != (QELT) 0 )
  1748.     {
  1749.     if( (k & 0x7fff) == (QELT) 0x4000 )
  1750.         {
  1751.         /* check all less significant bits */
  1752. #if WORDSIZE == 32
  1753. /* 0 1 2  3   4   5   6 */
  1754. /* S E M 1,2 3,4 5,6 7,8 */
  1755.         for( j=M+5; j<=NQ; j++ )
  1756. #else
  1757. /* 0 1 2 3 4 5 6 7 8 9 10 */
  1758. /* S E M 1 2 3 4 5 6 7 8 */
  1759.         for( j=M+9; j<=NQ; j++ )
  1760. #endif
  1761.             {
  1762.             if( ac1[j] != 0 )
  1763.                 goto yesrnd;
  1764.             }
  1765. /* round to even */
  1766.         if( (k & 0x8000) == (QELT) 0 )
  1767.           goto nornd;
  1768.         }
  1769. yesrnd:
  1770.     qclear( ac2 );
  1771.     ac2[NQ] = 0;
  1772. #if WORDSIZE == 32
  1773.     ac2[M+4] = 0x8000;
  1774. #else
  1775.     ac2[M+8] = 0x8000;
  1776. #endif
  1777.     addm( ac2, ac1 );
  1778.     if( ac1[2] )
  1779.         {
  1780.         shdn1(ac1);
  1781.         i += 1;
  1782.         }
  1783.     if( (i == 0) && (ac1[3] & 0x8000) )
  1784.         i += 1;
  1785.     }
  1786.  
  1787. nornd:
  1788.  
  1789. SC = 1;
  1790. shift( ac1 );
  1791. *e |= i;    /* high order output already has sign bit set */
  1792. p = &ac1[M+1];
  1793. #if WORDSIZE == 32
  1794. #ifdef MIEEE
  1795. ++e;
  1796. *e++ = *p >> 16;
  1797. *e++ = *p++;
  1798. *e++ = *p >> 16;
  1799. *e++ = *p++;
  1800. *e++ = *p >> 16;
  1801. *e++ = *p++;
  1802. *e++ = *p >> 16;
  1803. #endif
  1804. #ifdef IBMPC
  1805. *(--e) = *p >> 16;
  1806. *(--e) = *p++;
  1807. *(--e) = *p >> 16;
  1808. *(--e) = *p++;
  1809. *(--e) = *p >> 16;
  1810. *(--e) = *p++;
  1811. *(--e) = *p >> 16;
  1812. #endif
  1813. #else
  1814. /* 16-bit words */
  1815. #ifdef MIEEE
  1816. ++e;
  1817. *e++ = *p++;
  1818. *e++ = *p++;
  1819. *e++ = *p++;
  1820. *e++ = *p++;
  1821. *e++ = *p++;
  1822. *e++ = *p++;
  1823. *e++ = *p++;
  1824. #endif
  1825. #ifdef IBMPC
  1826. *(--e) = *p++;
  1827. *(--e) = *p++;
  1828. *(--e) = *p++;
  1829. *(--e) = *p++;
  1830. *(--e) = *p++;
  1831. *(--e) = *p++;
  1832. *(--e) = *p;
  1833. #endif
  1834. #endif
  1835. return 0;
  1836. }
  1837.  
  1838.  
  1839.  
  1840.  
  1841. /*                        qtoasc.c    */
  1842. /* Convert q type number to ASCII string */
  1843.  
  1844. /* Get values for powers of ten.  */
  1845. #include "qtens.h"
  1846.  
  1847. int qtoasc( q, string, ndigs )
  1848. QELT q[];
  1849. char *string;
  1850. int ndigs;
  1851. {
  1852. long digit;
  1853. QELT x[NTT], xt[NTT];
  1854. QELT *p, *r, *ten, *tenth;
  1855. QELT sign;
  1856. int i, k, expon;
  1857. char *s, *ss;
  1858.  
  1859. qmov( q, x );
  1860. sign = x[SIGNWORD];
  1861. x[SIGNWORD] = 0;
  1862. expon = 0;
  1863. ten = &qtens[NTEN][0];
  1864.  
  1865. i = qcmp( qone, x );
  1866. if( i == 0 )
  1867.     goto isone;
  1868. if( x[1] == 0 )
  1869.     {
  1870.     qclear( x );
  1871.     goto isone;
  1872.     }
  1873.  
  1874. if( i < 0 )
  1875.     {
  1876.     k = MAXNTEN;
  1877.     p = &qtens[0][0];
  1878.     qmov( qone, ac4 );
  1879.     qmov( x, xt );
  1880.     while( qcmp( ten, x ) <= 0 )
  1881.         {
  1882.         if( qcmp( p, xt ) <= 0 )
  1883.             {
  1884.             qdiv( p, xt, xt );
  1885.             qmul( p, ac4, ac4 );
  1886.             expon += k;
  1887.             }
  1888.         k >>= 1;
  1889.         if( k == 0 )
  1890.             break;
  1891.         p += NTT;
  1892.         }
  1893.     qdiv( ac4, x, x );
  1894.     }
  1895. else
  1896.     {
  1897.     k = MINNTEN;
  1898.     p = &qmtens[0][0];
  1899.     r = &qtens[0][0];
  1900.     tenth = &qmtens[NTEN][0];
  1901.     while( qcmp( tenth, x ) > 0 )
  1902.         {
  1903.         if( qcmp( p, x ) >= 0 )
  1904.             {
  1905.             qmul( r, x, x );
  1906.             expon += k;
  1907.             }
  1908.         k /= 2;
  1909. /* Prevent infinite loop due to arithmetic error: */
  1910.         if( k == 0 )
  1911.             break;
  1912.         p += NTT;
  1913.         r += NTT;
  1914.         }
  1915.     qmuli( ten, x, x );
  1916.     expon -= 1;
  1917.     }
  1918.  
  1919. isone:
  1920. qifrac( x, &digit, x );
  1921. /* The following check handles numbers very close to 10**(2**n)
  1922.  * when there is a mistake due to arithmetic error.
  1923.  */
  1924. if( digit >= 10 )
  1925.     {
  1926.     qdiv( ten, x, x );
  1927.     expon += 1;
  1928.     digit = 1;
  1929.     }
  1930. s = string;
  1931. if( sign != 0 )
  1932.     *s++ = '-';
  1933. else
  1934.     *s++ = ' ';
  1935. *s++ = (char )digit + 060;
  1936. *s++ = '.';
  1937. if( ndigs < 0 )
  1938.     ndigs = 0;
  1939. if( ndigs > NDEC )
  1940.     ndigs = NDEC;
  1941. for( k=0; k<ndigs; k++ )
  1942.     {
  1943.     qmuli( ten, x, x );
  1944.     qifrac( x, &digit, x );
  1945.     *s++ = (char )digit + 060;
  1946.     }
  1947.  
  1948. *s = '\0';    /* mark end of string */
  1949. ss = s;
  1950.  
  1951. /* round off the ASCII string */
  1952.  
  1953. qmuli( ten, x, x );
  1954. qifrac( x, &digit, x );
  1955. if( digit > 4 )
  1956.     {
  1957. /* Check for critical rounding case */
  1958.     if( digit == 5 )
  1959.         {
  1960.         if( qcmp( x, qzero ) != 0 )
  1961.             goto roun;    /* round to nearest */
  1962.         if( (*(s-1) & 1) == 0 )
  1963.             goto doexp;    /* round to even */
  1964.         }
  1965. roun:
  1966.     --s;
  1967.     k = *s & 0x7f;
  1968. /* Carry out to most significant digit? */
  1969.     if( k == '.' )
  1970.         {
  1971.         --s;
  1972.         k = *s & 0x7f;
  1973.         k += 1;
  1974.         *s = k;
  1975. /* Most significant digit rounds to 10? */
  1976.         if( k > '9' )
  1977.             {
  1978.             *s = '1';
  1979.             expon += 1;
  1980.             }
  1981.         goto doexp;
  1982.         }
  1983. /* Round up and carry out from less significant digits. */
  1984.     k += 1;
  1985.     *s = k;
  1986.     if( k > '9' )
  1987.         {
  1988.         *s = '0';
  1989.         goto roun;
  1990.         }
  1991.     }
  1992.  
  1993. doexp:
  1994.  
  1995. sprintf( ss, "E%d", expon );
  1996. return 0;
  1997. }
  1998.  
  1999.  
  2000. /*  QELT a[NQ], b[NQ];
  2001.  * qcmp( a, b )
  2002.  *
  2003.  *  returns +1 if a > b
  2004.  *         0 if a == b
  2005.  *        -1 if a < b
  2006.  */
  2007.  
  2008. int qcmp( p, q )
  2009. register QELT *p, *q;
  2010. {
  2011. QELT r[NQ];
  2012. register int i;
  2013. int msign;
  2014.  
  2015. if( ( p[E] <= (QELT) NBITS)  && ( q[E] <= (QELT) NBITS ) )
  2016.     {
  2017.     qsub( q, p, r );
  2018.     if( r[E] == 0 )
  2019.         return( 0 );
  2020.     if( r[SIGNWORD] == 0 )
  2021.         return( 1 );
  2022.     return( -1 );
  2023.     }
  2024.  
  2025. if( p[SIGNWORD] != q[SIGNWORD] )
  2026.     { /* the signs are different */
  2027.     if( p[SIGNWORD] == 0 )
  2028.         return( 1 );
  2029.     else
  2030.         return( -1 );
  2031.     }
  2032.  
  2033. /* both are the same sign */
  2034. if( *p == 0 )
  2035.     msign = 1;
  2036. else
  2037.     msign = -1;
  2038.  
  2039. i = NQ;
  2040. do
  2041.     {
  2042.     if( *p++ != *q++ )
  2043.         {
  2044.         goto diff;
  2045.         }
  2046.     }
  2047. while( --i > 0 );
  2048.  
  2049. return(0);    /* equality */
  2050.  
  2051.  
  2052.  
  2053. diff:
  2054.  
  2055. if( (unsigned int) *(--p) > (unsigned int) *(--q) )
  2056.     return( msign );        /* p is bigger */
  2057. else
  2058.     return( -msign );    /* p is littler */
  2059. }
  2060.  
  2061.  
  2062.  
  2063. /*
  2064. ;                                ASCTOQ
  2065. ;        ASCTOQ.MAC        LATEST REV: 11 JAN 84
  2066. ;                    SLM, 3 JAN 78
  2067. ;
  2068. ;    Convert ASCII string to quadruple precision floating point
  2069. ;
  2070. ;        Numeric input is free field decimal number
  2071. ;        with max of 15 digits with or without
  2072. ;        decimal point entered as ASCII from teletype.
  2073. ;    Entering E after the number followed by a second
  2074. ;    number causes the second number to be interpreted
  2075. ;    as a power of 10 to be multiplied by the first number
  2076. ;    (i.e., "scientific" notation).
  2077. ;
  2078. ;    Usage:
  2079. ;        asctoq( string, q );
  2080. */
  2081.  
  2082. int asctoq( s, y )
  2083. char *s;
  2084. QELT *y;
  2085. {
  2086. QELT yy[NQ+1], qt[NQ];
  2087. int esign, nsign, decflg, sgnflg, nexp, exp, prec;
  2088. QELT *p;
  2089.  
  2090.  
  2091. nsign = 0;
  2092. decflg = 0;
  2093. sgnflg = 0;
  2094. nexp = 0;
  2095. exp = 0;
  2096. prec = 0;
  2097. qclear( yy );
  2098. yy[NQ] = 0;
  2099.  
  2100. nxtcom:
  2101. if( (*s >= '0') && (*s <= '9') )
  2102. {
  2103. if( (prec == 0) && (decflg == 0) && (*s == '0') )
  2104.     goto donchr;
  2105. if( prec < NDEC )
  2106.     {
  2107.     if( decflg )
  2108.         nexp += 1;    /* count digits after decimal point */
  2109.     shup1( yy );    /* multiply current number by 10 */
  2110.     qmovz( yy, ac2 );
  2111.     shup1( ac2 );
  2112.     shup1( ac2 );
  2113.     addm( ac2, yy );
  2114.     qclear( ac2 );
  2115.     ac2[OMG+1] = *s - '0';
  2116.     addm( ac2, yy );
  2117.     }
  2118. prec += 1;
  2119. goto donchr;
  2120. }
  2121.  
  2122. switch( *s )
  2123. {
  2124. case ' ':
  2125.     break;
  2126. case 'E':
  2127. case 'e':
  2128.     goto expnt;
  2129. case '.':    /* decimal point */
  2130.     if( decflg )
  2131.         goto error;
  2132.     ++decflg;
  2133.     break;
  2134. case '-':
  2135.     nsign = -1;
  2136. case '+':
  2137.     if( sgnflg )
  2138.         goto error;
  2139.     ++sgnflg;
  2140.     break;
  2141. case '\0':
  2142. /* For Microware OS-9 operating system: */
  2143. #ifndef OSK
  2144. case '\n':
  2145. #endif
  2146. case '\r':
  2147.     goto daldone;
  2148. default:
  2149. error:
  2150.     printf( "asctoq conversion error\n" );
  2151.     qclear(y);
  2152.     return 0;
  2153. }
  2154. donchr:
  2155. ++s;
  2156. goto nxtcom;
  2157.  
  2158.  
  2159. /*        EXPONENT INTERPRETATION */
  2160. expnt:
  2161.  
  2162. esign = 1;
  2163. exp = 0;
  2164. ++s;
  2165. /* check for + or - */
  2166. if( *s == '-' )
  2167.     {
  2168.     esign = -1;
  2169.     ++s;
  2170.     }
  2171. if( *s == '+' )
  2172.     ++s;
  2173. while( (*s >= '0') && (*s <= '9') )
  2174.     {
  2175.     exp *= 10;
  2176.     exp += *s++ - '0';
  2177.     }
  2178. if( esign < 0 )
  2179.     exp = -exp;
  2180.  
  2181. daldone:
  2182. nexp = exp - nexp;
  2183.  
  2184. if( normlz(yy) )
  2185.     {
  2186.     qclear(y);
  2187.     return 0;
  2188.     }
  2189.  
  2190. yy[E] = EXPONE - 1 + NBITS - SC;
  2191. qmov( yy, y );
  2192. y[SIGNWORD] = nsign;
  2193.  
  2194. /* multiply or divide by 10**NEXP */
  2195. if( nexp == 0 )
  2196.     goto aexit;
  2197. esign = 0;
  2198. if( nexp < 0 )
  2199.     {
  2200.     esign = -1;
  2201.     nexp = -nexp;
  2202.     }
  2203.  
  2204. p = &qtens[NTEN][0];
  2205. exp = 1;
  2206. qmov( qone, qt );
  2207.  
  2208. do
  2209.     {
  2210.     if( exp & nexp )
  2211.         qmul( p, qt, qt );
  2212.     exp <<= 1;
  2213.     p -= NQ;
  2214.     }
  2215. while( exp < 8192 );
  2216.  
  2217. if( esign < 0 )
  2218.     qdiv( qt, y, y );
  2219. else
  2220.     qmul( qt, y, y );
  2221. aexit:
  2222. return 0;
  2223. }
  2224.